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Abstract. Given n points in Euclidean space E d , we propose an algebraic algorithm to 
compute the best fitting (d — l)-cylinder. This algorithm computes the unknown direction 
of the axis of the cylinder. The location of the axis and the radius of the cylinder are 
deduced analytically from this direction. Special attention is paid to the case d = 3 when 
n = 4 and n = 5. For the former, the minimal radius enclosing cylinder is computed 
■ algebrically from constrained minimization of a quartic form of the unknown direction of 

the axis. For the latter, an analytical condition of existence of the circumscribed cylinder is 
given, and the algorithm reduces to find the zeroes of an one unknown polynomial of degree 
at most 6. In both cases, the other parameters of the cylinder are deduced analytically. 
The minimal radius enclosing cylinder is computed analytically for the regular tetrahedron 
and for a trigonal bipyramids family with a symmetry axis of order 3. 

> 

in . 

£\\ . Keywords: Best fitting cylinder; smallest enclosing cylinder; minimal cylinder; circum 

\ scribed cylinder through five points; numerical algorithm 

00 

o 
o . 

2010 MSC codes: 51M04, 51N15, 65D10, 65K05, 90C26 



X! 

j_J ■ 1. Introduction 



Let E d be the d-dimensional Euclidean space, XJ d ~i be a (d — j)-dimensional affine 
subspace, and p a non negative real number. The set of points x of E d lying at distance 
smaller or equal to p from their orthogonal projection on U i is called a j-cylinder. We 
consider n data points in E d . A j-cylinder C containing n points is called a circumscribed 
j-cylinder if all n points are lying on the boundary of C. We consider further only (d— 1)- 
cylinders in this paper. A (d — l)-cylinder is defined with an axis, i.e. a unit vector u 
defining the direction of the axis and a point c of E d locating the axis, and a radius p. 

Given n data points Xi,i = 1, ...,n in E d , several cylinder problems are defined. 

Best fitting cylinder problem: Find the id— l)-cylinder (u, c, p) minimizing the standard 
deviation of the population of the n squared distances between the Xi and their projections 
on the axis, p 2 being the mean of this population. If it happens that the minimal standard 
deviation is null, all distances are equal to p and the resulting cylinder is circumscribed 
to the n points. 

Smallest enclosing cylinder problem: Find the (d— l)-cylinder (u, c, p) of minimal radius 
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p enclosing the n points. 

Circumscribed cylinder problem: Find the (d— l)-cylinder (tt, c, p) circumscribed to the 
n points, and if it is not unique, exhibit the one of minimal radius. The circumscribed 
cylinder may not exist. 

The circumscribed cylinder problem in E 3 seems to have appeared in the literature in 
1977 [1], but practical algorithms appeared much later [4], at the occasion of the analysis of 
circular cylinders through four or five points. Unfortunately, the practical computation of 
solutions is not trivial and requires the use of specialized packages. Recently, a numerical 
solver based on second order cone programming methods has been devised to compute 
smallest enclosing cylinders [15]. Bounds for the radius of minimal enclosing (i-simplices 
have been given [3], and the complexity of the smallest enclosing cylinder problem has 
been investigated for metrology applications [13]. 

In this paper we propose simple algebraic algorithms to solve the three cylinders prob- 
lems, which can be implemented without the use of sophisiticated packages. The moti- 
vation comes from recent molecular modeling studies in which the shape of many small 
molecules is expected to be suitably fitted with cylinders [9]. In this context, d = 3, and 
the smallest enclosing cylinder problem is the basic one we need to solve. We face to small 
n values, so that we are interested by the algebraic aspects rather than by the compu- 
tational complexity. Solving the best fitting cylinder problem is useful to introduce our 
approach, and solving the circumscribed cylinder problem is part of our full algorithm. 

2. The best fitting cylinder problem 

The usual scalar product of x and y is x'y, the quote denoting a vector of matrix 
transposition. The norm of x is ||x|| = {x'x) 1 / 2 . We define X as the rectangular array of n 
lines and d columns containing x\ at the line i, i = 1, ...,n. The squared distance between 
Xi and its projection on the axis is Aj = (xj — c)'{xi — c) — ((xi — c)'u) 2 . The mean of the 
population of the Aj is A, and the variance to be minimized is ^ 2~^(^« — A) 2 . 

Any point c lying on the axis of the cylinder can be used. Thus we decide to retain c 
as the projection of the origin on the axis, i.e. c'u = 0. We are left with the minimisation 
of a polynomial of u and c, subject to the orthogonality condition c'u = and to the 
normalization u'u = 1. Using the Lagrange multipliers K and L, the function to be 
minimized is 



For clarity we assume further than the mean of the n points is translated to the origin. 
Setting Vi = XiX^, the inertia matrix is T = ^ Vi and the covariance matrix is V = T/n. 
We also further assume that the set of the n points is not subdimensional, i.e. X is of 
full rank d, and T = X'X and V are invertible. The identity matrix of rank d is I, we 
define B 4 = I ■ Tr(Vi - V) - (V - V). Then, A, = Tr(Vi) - 2c'x; + c'c - u'ViU, and the 
expression of F cu expands as 




(1) 



h XX — A) 2 — Kdu — L{u'u — 1) 



\ J2(u'B iU - 2c'xi) 2 - Kc'u - L(u'u - 1) 



(2) 



The components of u and c are the 2d unknowns in (2). The unknown radius is p 2 = A. 
It is computed from c and u. The gradient of F cu relative to c is 
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G, 



,'Biu) + 8Vc - Ku 



(3) 



The stationary point is such that G c = 0, so that 



8c = V- 1 [Ku + ^J2(xiu'B i u)} 



(3) 



The Lagrangian K is expressed from setting u'c = in (3). Then, the analytical 
expression of the optimal c value can be used in (2), so that the minimization in (1) 
reduces to a non convex optimization problem in which there are only d unknowns, i.e. 
the components of u. Standard algebraic solvers can be used [7]. Moreover, starting points 
are random unit vectors u, which in fact are d—1 independant parameters values. This 
is beneficial in term of global optimization cost. E.g., for d = 3, a reasonable amount of 
random initial unit vectors should suffice to locate the global minimum of F cu , and there 
is no need of a huge of initial 6-tuples. 

The case where the solution is F* u = is of special interest. In this situation, all the 
n points are lying at distance p from the axis, and we face to a circumscribed cylinder 
problem. We have n equations Aj = A, i.e. 



Only Ti—l are independant equations because ^ Bi = and Yl x % = 0- Adding to them 
the two equations c'u = and u'u = 1, we get a system of n + 1 independant equations 
of 2d unknowns c and u. This system is determined when n = 2d — 1. When d = 3, we 
retrieve that at most a finite number of cylinders is expected to pass through 5 points in 
general position. More will be said further. 

3. The smallest enclosing cylinder problem 

A cylinder is a convex set, so that the smallest cylinder enclosing n points is the one 
which encloses the vertices of the convex hull of the n points. Thus, deleting the points 
interior to the hull before computing the smallest enclosing cylinder can drastically reduce 
the computational cost. Convex hull algorithms in 2, 3, and more dimensions have been 
described [5, 12]. 

When d = 3, the smallest enclosing cylinder problem is here the minimal radii enclosing 
cylinder problem, which should not be confused with the minimal height enclosing cylinder 
problem. This latter is a special case of the smallest enclosing 1-cylinder problem although 
the former is a smallest enclosing (d — l)-cylinder problem. The minimal height enclosing 
cylinder is defined by the closest enclosing parallel planes, and we know from proposition 
3 in [3] that these latter planes contain together at least one 4-tuple of vertices of the 
convex hull. The radius of the minimal height enclosing cylinder is then the radius of 
the smallest circle containing the projections of the points on one of the closest enclosing 
parallel planes. 

We return back to the minimal radii enclosing cylinder problem. Owing to the result 
given at the end of section 2, the smallest enclosing cylinder must be seeked in the set 
of minimal radius cylinders circumscribed to k points, k taking successively increasing 
values from 1 to 2d — 1. The smallest k value for which there is at least one minimal 
radius circumscribed cylinder enclosing all n points is retained. If there are several such 
circumscribed cylinders, the smallest radius one is retained. It may be not unique. 



vl 'BiU = 2c' Xi 



(4) 
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When d = 3, we look for minimal radius cylinders circumscribed to successively 1, 2, 

3, 4 and 5 points. Cases k = 1 and k = 2 are degenerate. The case k = 3 is in fact 
subdimensional: find the two closest lines enclosing 3 points. Still from proposition 3 in 
[3] we deduce that the associated minimal radius is half of the smallest of the three heights 
of the triangle defined by the three points. Only the cases k = 4 and k = 5 are non trivial. 
We consider them hereafter, assuming d = 3. 

4. The smallest cylinder passing through 4 points 

Here we set n = 4. For convenience, we set 7 = 2c and we define the vector b having 
four components 6j = u'BiU. The four equations in (4) are rewritten in matricial form 
X'y = b, and because we already assumed X to be of full rank, we express 7 from u in (5). 

7 = T- x X'b (5) 

The squared radius p 2 = A = Tr{V) — u'Vu + c'c is to be minimized under the 
constraints u'c = and u'u = 1. Setting W = I ■ Tr(T) — T and using the Lagrangians K 
and L, the objective function of u to be minimized is 

F = u'u ■ u'Wu + 7 ; 7 - 2K-i'u + L(l - u'u) (6) 

The case of the regular tetrahedron is solved analytically. 

Theorem 1. (a) There are three absolute minimal radius circumscribed cylinders to 
the regular tetrahedron. The ratio of the minimal radius to the length of the edge is 1/2. 
The axis of the three cylinders intersect at the center of the tetrahedron and are mutually 
orthogonal, following the directions parallel to the edges of the cube having for vertices the 
four ones of the tetrahedron plus the four ones got by reflection of these latter through the 
center, (b) There are six absolute maximal radius circumscribed cylinders to the regular 
tetrahedron. The ratio of the maximal radius to the length of the edge is 3\/2/8. The axis 
of the six cylinders following the directions parallel to the edges of the tetrahedron. 

Proof. We use: x' x = (1, 1, 1), x' 2 = (1, -1, -1), x' 3 = (1, -1, 1), x' 4 = (-1, -1, 1). 

Part (a): the center is at the origin and the covariance matrix V is the identity matrix 
/. The squared radius to be minimized is p 2 = Tr{V) — u'Vu + c'c = 2 + c'c, under the 
constraints u'c = and u'u = 1, with 7 = 2c = T~ l X'b. The absolute minimum at c = 
and p 2 = 2 will be valid under the condition that the four equations u'B^u = ^'xi = 
are indeed satisfied for some unit vector u. Then, u'B{U = Tr(Vi ~ V) ~ u'(Vi — V)u, i.e. 
u'BiU = x[xi — 3 — (u'xi) 2 + 1 = 1 — (u'xi) 2 . The solutions u are such that u'xi = 1 or 
u'xi = —1, which is written in matricial form Xu = 5, where 5 is a vector with each of 
its four component taking independantly the value 1 or the value —1, and 8'S = 4. So, 
u = T~ 1 X'5 = X'S/4 and u'u = 1. We deduce that the set of solution unit vectors u 
reduces to the set of the canonical base vectors and their opposite. The rest of the proof 
of part (a) is trivial. 

Part (b): the squared radius to be maximized is still p 2 = 2 + c'c. Noticing that XX' /A 
is the centering operator / — 11' /n, where / is now the identity matrix of size n = 4 and 
1 is the vector having its four components equal to 1, c'c = b' XT~ 2 X'bjA = 6'6/16. From 
(4), and because u'B^u = 1 — (u'xi) 2 , we have b'b = ^[1 — (u'xi) 2 ] 2 = [Yl{ u ' x i) 4 ] ~ 4. 
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After expansion: dc = [1 — (uf + u\ + ii|)]/2. The orthogonality constraint reduces to 
^(n'xj) 3 = 2AU1U2U3 = 0. Either ui or U2 or U3 is null. We solve the optimization 
problem in respect to the unkown vector (ttf , ^3)') which describes the boundary of the 
equilateral triangle having the three canonical base vectors as extreme points. The norm 
of this latter unknown vector is minimized when its extremity lies on the mid of any of the 
side of the equilateral triangle: uf = 1/2, = 1/2, ti| = 0, or u\ = 1/2, n| = 0,ii§ = 1/2, 
or u\ = 0,1*2 = 1/2, u\ = 1/2. Thus the unit solutions vectors u are in the directions 
of the edges and the maximal radius is 3/2. Note that minimizing c'c let to retrieve the 
canonical base vectors as solutions, reproving part (a). 

The results of theorem 1 are in agreement with those in [3, 4, 8] and with those for 
equifacial tetrahedra [2, 14], although the proof given here differs from the previously 
published ones. 

Returning to the case of the general tetrahedron, it is pointed out that we are seek- 
ing for an unknown direction rather than for an unit vector, so that the normalization 
condition can be weakened. The leading term u'u in (6) permits to write F as an ho- 
mogeneous quartic to be minimized under the orthogonality constraint j'u = 0. When 
this orthogonality constraint is satisfied by some vector u, it is still satisfied by any vector 
colinear to u, and when the gradient of the homogeneous quartic vanishes for some vector 
u, it vanishes again for any vector colinear to u. So, any non null vector u got during the 
Newton iterations can be renormalized without affecting the convergence of the iterative 
process. Such technique was successfully used in a minimal variance computation encoun- 
tered in the geometric docking problem [10, 11]. Some other practical implementation 
details follow. Random starting unit vectors are generated from an isotropic distribution. 
At the iteration k, the Lagrangians values are computed to minimize the norm of the full 
gradient of F. Then the Newton step sj, = —akH k l gj. is computed, where and Hy. are 
respectively the gradient of F at u k and its Hessian at Uk , and is a positive real, which 
is usually set to 1 when the quadratic convergence is achieved (see also the discussion at 
the end of the paper). However the resulting vector + does not satisfy anymore to 
the constraints. The orthogonality condition can be restored via solving the minimization 
problem 

Miri{ s y(s - s k )'(s - s fc ), subject to g' k s < 77, and to u k+ i^k+i = (7) 

The parameter rj in the first constraint is set to ensure that, using a first order approx- 
imation, there is a sufficient decrease of F. For that, we should impose a negative value 
to g' k s, and we simply set rj = g' k Sk- But solving (7) needs a costly extra work due to the 
orthogonality constraint. Thus, we expand the expression of the orthogonality constraint 
and we approximate it to the first order of the unknown s. The modified constrained 
minimization problem (7) is solvable analytically, and the Newton step s is computed at 
low cost. Then, as mentioned above, the normalization constraint is restored via setting 
\\ u k+l\\ = 1- 

Our implementation of the solver was programmed in f77 with double precision (i.e. 
64 bits) floating numbers. The termination criterion was \cos(u, c)| < 10~ 10 and |(i*fc+i — 
F; c )\/(Fk+i + Fk) < 10 -10 . Alternate termination criterions may consider the norm of the 
gradient, the norm of the Newton step s, etc. 

Thousands of random tetrahedra were generated from various probability laws, includ- 
ing quasi-flat and quasi-linear sets. The implementation of the Newton method appeared 
to be effective for most random initial unit vectors. In fact, for most tetrahedra, the 
convergence was observed for all initial vectors. In the worst situations, reducing signifi- 
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cantly the length of the Newton step at each iteration before restoring the orthogonality 
constraint and operating with 100-250 initial values suffice to overcome the problem of 
convergence failure and to observe a significant number of times the convergence. Any- 
way, there may be several locally minimal radius cylinders enclosing a tetrahedra, and the 
upper bound of the number of local minima has been shown to be 9 [4]. So, we indeed 
needed to operate with a sufficient number of random initial unit vectors. 

A typical example is the regular tetrahedron x' x = (1,1,1), x' 2 = (1,-1,-1), x' 3 = 
(1, —1, 1), x' A = (—1, —1, 1). Two locally minimal radii were found: p = 1.41421356 and 
p = 1.50000000. Among the 100 performed minimizations, the smallest radius was found 
32 times and the other one was found 68 times. The numbers of iterations ranged from 4 
to 35, with a mean value around 9.7. The axis of the cylinder associated to the smallest 
radius is in agreement with the one given in theorem 1. We observed that the other radius 
corresponds to the six unit vectors in the six respective directions defined by the edges of 
the tetrahedron, for which the calculated value is indeed p = 3/2. 



5. The circumscribed cylinder problem 

Here we set n = 5. The convex hull of the 5 points is either a tetrahedron with one of 
the 5 points lying in its interior, or a trigonal bipyramid. Falling in this latter configuration 
is a necessary condition of existence of the circumscribed cylinder, but it does not suffice. 
According to (4), we must find the zeroes of the system (8) of 6 unknowns c and u 

u'BiU = 2c' Xi ,i = l, 5; u'c = 0; u'u = 1 (8) 

The system (8) contains only 6 independant equations because the 5 first ones are 
dependant (sum to zero). Again, the radius is computed from c and u and the center 
c = 7/2 is computed according to (5): 

p 2 = A = Tr{V) - u'Vu + c'c 
7 = T~ l X'b, hi = u'B iU 

Thus, reporting the expression of c in some linear combination t of the 5 first equations 
of (8) lead to a system of only 3 equations of 3 unknowns, i.e. the components of u. We 
define the vector 1 as a vector having its 5 components equal to 1. A being centered and 
of full rank, there is only one direction orthogonal to the 4-dimensional subspace defined 
by the the three independant columns of X and the vector 1, this latter being itself 
orthogonal to each of the 3 columns of X, i.e. X'l = 0. We set the linear combination t 
to follow this unique direction, and we conventionnally normalize t't = n, and a plus sign 
is arbitrarily attributed to the component of t having the largest absolute value. In other 
words, the computation of t needs that we perform a Gram-Schmidt orthogonalization of 
the columns of X and of the vector 1. Alternatively, an adequate rotation can be made 
so that X is positioned in its principal components set of axis (i.e. T = X'X is diagonal), 
so that the columns of X and the vector 1 are all orthogonal, thus providing an easy way 
to compute t. 

i=5 i=5 

Having X't = means that ^ Xj = 0. Then we define the matrix M = ^ UBi, and 

i=l i=l 

we consider the system (9) issued from the linear combination t of the 5 first equations of 
(8): 
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u'Mu = 0; u'c = 0; u'u = 1 (9) 

This system, where c is a quadratic function of u, has three equations of three unknowns. 
It is valid if and only if M is not null. 

Theorem 2. X being full dimensional, the matrix M is null if and only if two points 
are identical. 

Proof We assume M = £ fcflj = 0. We have t'l = 0, and so, £ • Tr(V;) - V^) = 0. 
Because Tr(M) = 0, we have also ^UTriyi) = and thus J2^i x i x 'i = 0- We define the 
diagonal matrix G with its diagonal element ©jj = ti, i = 1,...,5. Then, X'QX = 0, so 
that we can write the matrix product 



X' 



X'Q' 



X QX 



T 








x'e' ex 



The result of the matrix product above is a square matrix which cannot be of full rank 
because it is the product of a 6 columns and 5 lines matrix by its transposed. Because 
X is of full rank, the block X'Q'QX = ^ tfxixi' must have at least one null eigenvalue. 
So, there is some vector uj such that uj' X'Q'QXuj = 0, and so QXuj = 0. We set X1234 
as the four lines array of the points x\, X2, X3, X4 , and 61234 as the associated diagonal 
submatrix of 0, and we remark that X1234 is of full rank because if it was not, there 
would be some vector orthogonal to x%, X2, £3 and X4, and thus it would be orthogonal 

i=4 

to £5 because X5 = — ^2 X{, thus implying that X is not of full rank, a contradiction. 

i=i 

This remark stands for all the submatrices of X containing four of the five points. Now, 
if ©1234 would be invertible, we would have found a vector uj such that X1234W = 0, which 
is impossible, so 61234 is not invertible, and there at least one zero element among t±, 
£3, £4. This conclusion is valid for all four-tuples of {1, 2, 3, 4, 5}, so that there are at least 
two distinct null elements among t±, t2, £3, t±,t§. 

Conventionnally, we set £4 = and t$ = without loss of generality. From the 
orthogonality conditions X't = and l't = 0, we have t\Xi + £2^2 + £32:3 = and 
^1 + ^2 + ^3 = 0, which means that xi, X2, x$ are aligned. Assuming t^t^ 7^ 0, 
X3 = (tiXi + t2X2)/(ti+t2) and expanding the expression of M = tiXix'i + t2X2x' 2 + £3X3X3 
leads to M = £1*2(^1 — %2)(x\ — X2) 1 /(h + £2)- The assumption M = would imply 
here that x\ = X2 = X3, which is impossible due to the full dimensionnality of X. So, 
ii^2^3 = 0, and we must assume that a third element of t is null. We set conventionnally 
£3 = 0, and we have both t\ + £2 = and t\X\ + £2^2 = 0, meaning that x\ = X2 because 
t is not null by definition, and M is indeed null. 

Conversely, we consider a full rank set such that two points, say 1 and 2, are identical. 
The unique direction t is colinear to (1, —1, 0, 0, 0) and thus M = 0. 

Having two identical points means that in fact we face to the four points problem, and 
this situation was considered in section 4. It may also be considered that suppressing the 
equation u'M u = lead to an underdetermined system, and so the problem of finding the 
zeroes of a function should be transformed in an optimization problem where the radius 
of the cylinder is to be minimized. Anyway, we further assume that we have never two 
identical points. 



Theorem 3. The set of solutions of (8) and the set of solutions of (9) are equal. 
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Proof. Obviously, the set of solutions of (8) is included in the set of solutions of (9). In 
order to ensure that the solutions of (9) satisfies to (8), we should check that they indeed 
satisfy to the five equations u'BiU = 2c' a^, i = 1, 5, or, in matricial form, b = 2Xc. We 
recall that the center is computed from 7 = 2c = T~ 1 X'b, with hi = u'BiU, i = 1, ...,5. 
The 5-dimensional vector b can be decomposed in the basis of 5 independant vectors build 
with the 3 columns of X and the vector 1 and the vector t, i.e. b = 2Xfix +/3il+/3ti, where 
fix contains the three coefficients associated to the columns of X and fix and fit are the 
coefficients associated respectively to 1 and t. The sum of the five Bi is null, thus l'b = 
and so fix = 0. From u'Mu = we have t'b = £ ti(u'B iU ) = 0. Thus, 2t'Xfi x + fitt't = 0, 
from which fit = because t'X = 0. So, b is indeed a linear combination of the columns of 
X, i.e. b = 2Xfix- The computed center is 2c = T~ l X'b, so fix = c and indeed b = 2Xc. 
Remark: the assumption M^O was not used in the proof above. 

Lemma 1. There are at most 6 cylinders circumscribed to 5 points. 

Proof. The homogeneous polynomial equations u'Mu = and u'c = in the system 
(9) define the intersection of two curves in the projective plane of respective degrees 2 and 
3. The normalization condition u'u = 1 is just a practical way to avoid a null vector u, 
and in fact all non null vectors colinear to a solution u at the intersection of u'Mu = and 
u'c = define the same intersection and lead to the same cylinder. So, applying Bezout 
theorem shows that we should have at most 6 cylinders passing through 5 points. This 
result was already proved in [4] via an other method. 

Lemma 2. There is no cylinder circumscribed to 5 points when either M or —M is 
positive definite. 

Proof. Obvious from u'Mu = 0. 

The matrix M = £ U(I • Tr{Vi) - Vj) = I-Tr(X'QX) -X'OX, characterizes a quadric 
surface. Except in the particular situation where an eigenvalue of M vanishes, the surface 
u'Mu = t'b = 0, i.e. ^t^u'u ■ x\xi — (u'xi) 2 ) = 0, is an elliptical cone. The cubic surface 
is it' 7 = 0, with 7 = T _1 ^2 Xiiu'u ■ x\xi — iu'xi) 2 ). 

Solving the system (9) can be reduced to an one unknown problem. We assume without 
loss of generality that the data set is rotated to have M being diagonal with eigenvalues 
Ml > > /^3- The case //1/X3 > was proved to correspond to no circumscribed cylinder 
and the case fi\ = ^3 = was proved to correspond to an infinite number of cylinders 
(four points problem). 

The equations u'Mu = and u'u = 1 define a linear system of two equations of the 
the three squared components of u. One of the components u\ or U2 or u$ is taken as 
the unknown parameter, and the expressions of the two other components are reported 
in n'7 = 0. Due to the free choice of the signs of these two other components, we have 
in fact four one unknown equations to solve. Assuming that we performed the rotation 
above and that the eigenvalues of M are separated, the relations between the components 
of u are 



u 2 2 = («i(M3 - Mi) - M3VO2 - A*a) «| = ( u i(^i ~ + M2VO2 - Ma) 
«i = («i(A*3 - M2) - M3)/ (pi ~ M3) uj = («!(Afc - Mi) + Mi)/(Mi - M3) 
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«i = («i(A*2 - Ma) - M2)/(Mi - M2) «1 = (^l(M3 - Mi) + Mi)/(Mi - Ma) 



The components of u must be in [0;1], but from the relations above we get tighter 
bounds for the unknown component to be selected. When [12 > 0: 

< u\ < 1 - Mi/(Mi - Ma) 
< «| < 1 - ^2/(M2 - M3) 
M2/ (M2 - M3) < u\ < 1 + ^3/ Oi - fl 3 ) 

And when [12 < 0: 

-M2/ (Ml - M2) < u\ < 1 - /xi/ (/xi - /z 3 ) 
< u 2 , < 1 + /x 2 /(mi - M2) 
< u 2 < 1 + ^ 3 /(Mi - /i 3 ) 

When either /ii = ^2 or ^2 = /x 3 , some of the relations and intervals above are invalid. 
The case fj,± = 113 is such than there is no cylinder to find, either because M or — M is 
positive definite, or because M = and we have a four points problem. Outside these 
situations, we can always retain U2 as the unknown parameter. Whatever unknown we 
select, either m, U2 or tt 3 , the resulting function is not a polynomial. 

Theorem 4. The system (9) can been solved via extracting the real roots of a polynomial 
of degree at most equal to 6. 

Proof. Having rotated the set so that M is diagonal, we rotate it again around the 
second eigenvector of M with an angle 02 to be specified. 

u'Mu = u\(iiicos 2 {ct2) + M3 s 2n 2 (a 2 )) + w 2 M2 + u 2 (/iisin 2 (a2) + ^3Cos 2 (a 2 )) 
-2niu 3 sin(a 2 )cos(a2)(Mi - M3) 

We consider the general case where < 0. The case where //1/U3 = and either 

[h\ 7^ or /i3 7^ will be considered later. 

We select 02 so that //isin 2 (a 2 ) + M3 cos2 ( a 2) = 0, i.e. tg(w) = a/— M3/Mi- The quadric 
becomes a parabola in the projective plane defined by the coordinate u\. 



Reporting u^/ui in the expression of the cubic curve in the projective plane leads to a 
polynomial of degree six in U2/U1, which has at most 6 real roots. We deduce u^/u\ from 
(11) and then we have the direction of u. So we get at most 6 cylinders. 

The existence of solutions such that u\ = should be checked. When ^2 7^ 0, we 
deduce from (10) that U2 = 0, and hence v! = (0,0, ±1) should satisfy to 1/7 = 0, so that 
the coefficient of m| in the cubic form is null. If the solution u' = (0, 0, ±1) is indeed found, 
the other ones found in the projective plane from (11) are such that the curve defined from 
1/7 = degenerates to a quintic of 112/ u\ and we cannot expect more than 5 real roots, 
i.e. there are at most 6 cylinders. 

When [12 = 0, the homogeneous cubic of U2 and U3 offers at most 3 real roots when 
solved either in U2/U3 or U3/112 (1*2 and U3 cannot be both null). The solutions other than 



u?(/ii + /U 3 ) + ul(i2 - 2u 1 u 3v /= MiM3 = 



(10) 
(11) 



( u 3\ — M2 ( U2 \2 _ /J1+M3 
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u' = (0,U2,us) found in the projective plane from (11) are such that the curve defined 
from it' 7 = is a cubic of uiju\ and we expect at most 3 real roots, i.e. there are at most 
6 cylinders. 

We consider now the case where /ii//3 = and either m ^ or ^3 7^ 0. Selecting 
«2 = so that M is diagonal, we have n\u\ + ^2^2 + /^3^| = 0. When [12 7^ 0, there is 
only one potential solution which is either v! = (0, 0, ±1) or v! = (±1, 0, 0). When ^2 = 0, 
either u\ = when fj,\ 7^ or u 3 = when 113 7^ 0. When fj,\ 7^ 0, it' 7 is an homogeneous 
cubic of it2 and 113 which offers at most 3 roots when solved either in ^2/^3 or U3/U2 («2 
and U3 cannot be both null), so that we expect at most 3 cylinders. A similar conclusion 
is got when /X3 / 0. 

In any case, the system (9) can be solved by finding the roots of a polynomial of degree 
at most 6, and we cannot get more than 6 cylinders, which redemonstrates lemma 1. 



6. An example of trigonal bipyramid with a symmetry axis of order 3 

We present the case of a bipyramid BV symmetric around its equilateral triangular basis 
with a symmetry axis of order 3 orthogonal to the triangular basis. We set x\ = (0, 0, h), 
x > 2 = (0,0, -h), 4 = (1,0,0), x' 4 = (-1/2,^3/2,0), 4 = (-1/2,-^3/2,0). The set is 
centered. 



/ ■ Tr(Vi) — V]_ = I • Tr(V 2 ) - V 2 



h 2 











h 2 















/ • Tr(V 3 ) - V 3 



I ■ Tr(V A ) - V A 



3/4 y/3/4 
V3/4 1/4 I I-Tr(V 5 )-V 5 



3/4 




(1/V6) 



/ 3 \ 

3 

-2 
-2 

v -2 y 



M = (1/V6) 



6h 2 - 3 
6h 2 - 3 
0-6 





1 

1 

-V3/4 

1/4 

1 



So, there is no circumscribed cylinder when h < \/2/2. We assume h > V2/2. 



Bi = Bo 



B 4 



Br 



(6/i 2 -3)/10 

(6/i 2 -3)/10 

-3/5 

(-4/i 2 -3)/10 

(7-4/i 2 )/10 

2/5 

(9-8/i 2 )/20 v/3/4 

y/3/4 (-l-8/i 2 )/20 

2/5 

(9-8/i 2 )/20 -\/3/4 

-V3/4 (-l-8/i 2 )/20 

2/5 
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b x = b 2 = 3(2/i 2 - 2 + u\ + u%)/10 
b 3 = (4 - u\{Ah 2 + 7) + u 2 (3 - 4/i 2 ))/10 

6 4 = (8 + (1 - 8/i 2 )u 2 - (9 + 8/i 2 )n| + 10V3um 2 ))/20 

6 5 = (8 + (1 - 8h 2 )ul - (9 + 8/i 2 )n| - 10\/3uiu 2 ))/20 

/ 2/3 -1/3 -1/3 \ 

T~ l X' = \/3/3 -V5/3 7 = 2c 

\ 1/2/i -l/2/i / 

Solving u'Mu=0 with 1/7 = 0, i.e. (2/i 2 + l)(uf + u 2 ,) = 2 and ui(uf - 3it|) = leads to 
the six desired directions ii and we easily deduce Theorem 5 below. 

Theorem 5. 

(a) The symmetric bipyramid BV of interapex half-height h has no circumscribed cylinder 
when h < v / 2/2. 

(b) When h > V2/2, BV has six circumscribed cylinders. The axis are in the directions 
u' = (0, ^2/(2h 2 + l),±.y/(2h 2 - l)/(2h 2 + 1)), which intersect aid = (l/(4/i 2 + 2), 0, 0), 
plus their rotated images of respective angles 2ir/3 and Att/2> around the interapex symmetry 
axis of order 3. All six cylinders have radius p = {Ah 2 + l)/(4/i 2 + 2). 

Neither the condition of existence of the cylinders circumscribed to BV nor their ana- 
lytical calculations seem to be previously reported in the literature. They are mainly 
consequences of Theorem 3 and Lemma 2, which themselves cannot be trivially deduced 
from existing results in cited papers and provide a simpler approach to the circumscribed 
cylinder problem in E s . 

7. Discussion and conclusion 

Some results in this work were already available in previous papers such as [2, 4, 14]. 
We had clearly mentioned these references at the appropriate places. Since these results 
are obtained here via different techniques and add both to the clarity and the self-content 
character of the paper, we left them. Having said that, we recall that several novelties are 
introduced in this paper. First, formulating the problems through variance minimizations 
led us in the three-dimensional case to the original systems (6) and (9), which offer partic- 
ularly simple expressions of the unknown axis u. It means that iterative numeric solvers 
need only random 3-tuples, and can be programmed by most users having some knowledge 
of constrained optimization, without needing to install specialized mathematical packages. 

In the case of the circumscribed cylinder to five points, we produced a simple and new 
condition of existence of the cylinder relying on the eigenvalues of the matrix M in section 
5. This matrix and its properties were yet unpublished. The explicit expressions of the 
polynomials of one variable at the end of section 5 are original, too. However, we found 
easier to programme directly the solver of (9) rather than seeking for the roots of the 
one variable polynomials. It was just an arbitrary developper's choice. Then, the third 
reviewer of this paper mentioned that the presence of square roots in the polynomials 
could be penalizing, and thus particular methods such as [6] could be used. 

Finally, the original analytical solution for the symmetric bipyramid in section 6, which 
relies on the explicit calculation of the matrix M, is outlined. 

By no way we claim that the original present approach is better than the previous 
ones, such as in [4], and the present software, named CYL, downloadable for free at 



-u( + u 2 )/2 

U\U2 
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http://petitjeanmichel.free.fr/itoweb.petitjean.freeware.html, is neither claimed to offer 
the best possible implementations of the methods nor is claimed to work better than 
other cylinder computation softwares. Moreover, the second reviewer informed us that the 
algebraic approach in [4] could guarantee the result, although our solver (see in particular 
section 4) cannot offer this guarantee for all datasets, despite the encouraging results we 
got. 

It would have been of high interest to compare our practical experiment results with 
those of existing methods such as the algebraic method in [4]. But, apart for the regular 
tetrahedron (detailed numerical results at the end of section 5), and for the bipyramid in 
section 6, we found only one public dataset, mentioned in [15]: we retrieved the optimal 
radius of the 12 points set minimal enclosing cylinder computed in [15] with all significant 
digits. 

We propose to those readers who have access to some cylinders computation software, 
to exchange with us any data and practical experimentation results they like. 
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